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Abstract 

In this work we characterize the configurational space of a short chain of colloidal particles as function of 
the range of directional and heterogeneous isotropic interactions. The individual particles forming the chain 
are colloids decorated with patches that act as interaction sites between them. We show, using computer 
simulations, that it is possible to sample the relative probability of occurrence of a structure with a sequence 
in the space of all possible realizations of the chain. The results presented here represent a first attempt to 
map the space of possible configurations that a chain of colloidal particles may adopt. Knowledge of such 
a space is crucial for a possible application of colloidal chains as models for designable self-assembling 
systems. 

PACS numbers: 64.70.pv, 64.75.Yz,64.60.De 
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I. INTRODUCTION 



Self-assembly is the process by which a substance exclusively driven by non-covalent inter- 
actions spontaneously develops into a specific long-lived conformation with a well defined struc- 
ture [1J. Self-assembling is common in natural bio-polymers, such as DNA, RNA and, in par- 
ticular, proteins, that exhibit exceptional self-assembling properties. Proteins, the fundamental 
building blocks of every living organisms, have the remarkable property that their structure and, 
therefore, their function is encoded in the one-dimensional sequence of the structural elements that 
compose them. Although proteins vary strongly in size, structure and function, they are all com- 
posed of the same 20 fundamental chemical units called amino acids. This protein alphabet insures 
an enormous variety of combinations, however only few sequences will have a well defined stable 
ground state, each of which consists of complex arrangements of predominantly three types of 
secondary structures: alpha helices, beta sheets and random coils (2]-|4]|. Hence, bio-polymers are 
a clear reference point for the development of artificial self-assembling systems based on modular 
subunits. 

Recently, it has been shown [5] that the minimal set of constraints imposed on configurational 
space by the hydrogen bonds along the protein backbone and the self-avoidance of the residues 
are sufficient to enable the successful folding of real protein structures. More generally, it is 
known that the specific geometry of the backbone and the directionality of the hydrogen bonds 
"pre-sculpt" the configurational space of proteins to the sub-space of currently known proteins 
structures [[6]-[9]], and allow for the design of sequences capable of folding back to their target 
structures 0. We believe that it is possible to extend the concept of configurational space pre- 
sculpting to generalized chains of particles capable of folding into geometries very different from 
those observed for natural proteins. However, for the future design of complex structures it is 
crucial to characterize the variety of configurations that a given set of interactions can generate. 
This is equivalent to the still ongoing process of mapping the "protein universe", where huge 
efforts are devoted to finding all protein structures corresponding to all sequences expressed in 
living organisms ifTOll . 

In this work, we explore the total phase space of colloidal patchy polymers. Such phase space 
is the ensemble of all chain configurations defined by a structure and a sequence along the chain 
of particles. Each particle is decorated by spots, or patches, with interactions that mimic the 
directionality of the hydrogen bonds. Such directional interactions confine the configurational 
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space of the synthetic protein analogue to a sub- space of compatible structures. We map the phase 
space for a range of values of the interaction potential and for two different forms of the interaction 
potential. We have chosen as subunits patchy particles, because in addition to mimicking the 
properties of hydrogen bonds, they have a rich ensemble of self-assembling properties lfTT1 - [T3ll . 
Experimentally, patchy colloidal particles are spherical particles with a hard core and a radius that 
varies from the nanometer to the micrometer scale and several methods have been developed for 
the experimental realization of such particles 03-16]. A beautiful review on the experimental 
methods to produce patchy particles can be found in the paper of Bianchi et al. lfT2l . 

The results presented in this paper represent a first attempt to link the configurational space of 
a generic chain of colloidal particles to simple geometrical parameters of the interaction potential 
between the particles. Such a space represents a "phase diagram 46 of the designability of the 
system, because the occurrence probability reflects the propensity of a generic sequence to adopt 
a given structure. Hence, the larger the space the more versatile is the corresponding model. 
Below, we will first introduce the details of the model and explain the computational methods 
used to sample the complex configurational space. Then, we will study the properties of chains of 
particles with three patches each, where each patch has different interaction ranges relative to the 
particle-particle interaction, as well as two different functional forms. 

II. METHODS 
A. Model 

The system we study here consists of a single chain of particles, each of which can occur in 
different types mimicking the various residues of a real protein. Each particle is decorated with 3 
patches, two of which are used to link the particle with neighboring particles and form the chain 
as illustrated in Fig. [TJ These connecting patches interact via a simple harmonic spring potential, 

£link(r) = V. (1) 

Here, r is the distance between the two anchoring patches and K = 5 k^TR^, where Rhc is the 
"hard core" radius of the colloids. In what follows all the energies are given in units of &#r re f, 
where T re f is a reference temperature that sets the scale of all interactions, hence in what follows 
the temperatures do not have units. 
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Figure 1: Real-space representation of the backbone of the patchy polymer. The patches are indicated by 
small white spheres, while the large turquoise spheres represent the hard core of the colloidal particles. If 
a comparison between this chain and the caterpillar backbone is made 0, it is easy to see how the patches 
are reminiscent of the O and H atoms of the protein backbone. The insets show a schematic representation 



The model colloidal polymer studied here is based on the caterpillar model [5j. Accordingly, 
the colloidal particles interact pairwise via a smoothed square well pair potential 



where the subscripts A and B refer to the type of two interacting particles and r is the distance 
between their centers. The parameter K = 0.2Rhc determines the width of the range over which 
the potential changes from ~ to ~ £ab, and r max is the distance at which £A#(/rmax) = £a#/2. 
We have varied r max in the range r max = 2.2,2.5,3.0,4.0,6.0 [Rhc]- The parameter £ab, which 
depends on the types A and B of the interacting particles, determines the depth of the well, i.e., 
the strength of the interactions. Depending on the particular type of the interacting particles, 
£ab can be positive or negative, leading to repulsive or attractive interactions, respectively. Each 
particle is then assigned a particular identity that will define the interactions with the other particles 
along the chain. There is no unique way to choose the values of the £ab, and we will show that 
with just an alphabet of 2 letters it is possible to produce a wide range of structures. The 2 
letters will only distinguish between hydrophilic (P) and hydrophobic particles (H) with £pp = 0, 



of the distance R between the patches and the alignment angles 0i and 02, which, according to Eq. [3] and 
Eq.|4], determine the interactions between patches not involved in the polymer backbone. 




(2) 
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Ehh = —MkBT YQ f, Zhp = — a/2^r re f. Here a is an adjustable parameter that we will fix to keep 
the second virial coefficient B^ B ~ 1 constant for the various ranges r max of the potential. The 
value of #2 B corresponds to the contribution of Eab in Eq.[2]to the total energy and was calculated 
with numerical integration for each value of r max . 

As patch-patch interaction, acting only among the free patches not involved in the links along 
the chain, we use two sets of directional interactions. The first one is inspired by the hydrogen 
bond interaction from the caterpillar model, which is represented by a 10-12 Lennard- Jones type 
potential with an additional angular dependence that takes into account the directionality of the 
patch bonds IH71 . 

12 /fTm\10l 



Ep\ = — 8pi (cos6i cos62) v 



(3) 



Here, v = 2 as given in ifTTL r is the distance between the patches, and 6i and 62 are angles 
between the patches and the centers of their corresponding spheres, as indicated in the left inset 
of Fig. [T] Op\ represents the position of the minimum of the potential; once a value for Gpi is 
fixed then also the range of the potential is determined. The smaller the value of Gpi the shorter 
the range of the potential will be. For each of the values of r max in Eq. [2] we considered short 
and long range directional interactions summarized below. As for a that defines the £ab scale 
factor in Eq. |2j we determined the parameter £pi by imposing that the second virial coefficient 
B% 1 of the directional potential in Eq. [3] is equal to the second virial coefficient B^ B for each pair 
(/max 5 <3pi). The hard core of the spheres ensures that only the maximum of angular term close to 
71 is accessible, i.e. —71 corresponds to configurations that are sterically inaccessible. 

The second type of directional potential is represented by the same square-well like potential 
defined in Eq. [2] modulated by the same angular dependence used for Ep\ in Eq. [3j but a different 
definition of the angles 61 and 62, 



Ep2 



-8p2 (cos 61 cos62) v 



1 



(c P2 -r)/K 



l+e 



if r > R HC , 

(4) 

tir<R HC9 



Here, v = 2, K = 0.2Rhc, r is now the distance between the centers of the spheres, Gp2 defines 
the range of the potential, and 61 and 62 are now the angles of the projection of each patches on 
the vector joining the centers of the interacting spheres, as indicated in the right inset of Fig. [TJ 
We choose a different definition for the distance r in the potential £>2, compared to Ep\, because 
we compute the same distance for the potential Eab i n Eq. [2j As the interaction P2 in Eq. [4] is 
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defined between the centers of the spheres, the range <Jp2 will be longer than Gp\ from Eq. [3] by 
exactly 2 colloids radius Rhc- Again, for each of the values of r max in Eq. [2] we considered short 
and long range directional interactions summarized in Table [I| The parameter £p2 is determined 
by imposing that the second virial coefficient B^ 2 of the directional potential in Eq. [4] is equal to 
the second virial coefficient B^ B for each pair (r max ,Gp2)- The directionality of the patch bonds, 
encoded in the angular term in Eqs.[3]and [4], are each essential to pre-sculpt the conformational 
space in analogy to the role of hydrogen bonds that characterize the secondary structure elements 
typical of proteins. 

We considered the directional interactions in Eqs.[3]and [4] as a way to estimate the dependence 
of configurational space on the different ways of imposing the directionality. Since the directional 
contribution as well as the window of interaction ranges considered are the same, the only differ- 
ence between the two interactions is in the shape of the potential. While the interaction in Eq. [3] 
has a sharp minimum, the square- well profile of the potential in Eq. |4} introduces less frustration 
when combined with the isotropic interaction of Eq.[2} which has the same functional form. More- 
over, there are also experimental reasons to consider the above mentioned interactions. The first 
potential can represent a short-range directional attraction on the surface (e.g. Hydrogen Bonds, 
short DNA), while the second is representative of a class of particles, where the patch is produced 
by creating a dent in the isotropic interaction. 



Parameter 


Setl 


Set 2 


Set 3 


Set 4 


Set 5 


Set 6 


Set 7 


Set 8 


?max [Rhc] 


2.2 


2.5 


3.0 


4.0 


4.0 


6.0 


6.0 


6.0 


a [k B T ie f\ 


2.6 


1.7 


1.0 


0.441 


0.441 


0.14 


0.14 


0.14 


Gpi [Rhc] 


0.2 


0.2 


0.2 


0.2 


1.0 


0.2 


1.0 


2.0 


£pi [A^ef] 


5.95 


5.95 


5.95 


5.95 


4.0 


5.95 


4.0 


3.1 


Gp2 [Rhc] 


2.2 


2.2 


2.2 


2.2 


3.0 


2.2 




4.0 


E[>2 [kaT^f] 


5.0 


5.0 


5.0 


5.0 


2.7 


5.0 




1.38 



Table I: Values of the potential parameters used in our simulations. It is important to remember that as the 
range C/>2 is calculated between the centers of the colloids it will be longer than o P \ by exactly 2 colloids 
radius Rhc 
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B. Simulation methodology 



The first step in the characterization of the patchy polymer model is the classification of the 
different structures according to how easy it is to design them. In order to do so, we will compare 
the designability of many chain configurations by estimating the number of sequences that fold into 
each structure. In order to explore this vast space we combine a particle identity mutation move 
(see below) with the standard set of pivot crankshaft and single particle moves used to explore the 
configurations of chains of particles lfT8l . For relatively short chains it is possible to identify the 
most recurrent structures and compare their designability using the algorithm described below. We 
considered first the simple case of a chain of 20 colloidal particles chosen from an alphabet of size 
2 as described in the section II.A. 

As we aim to imitate the designability property of natural proteins, we based the identity muta- 
tions move on the design scheme successfully used in protein studies [5j. As in the conventional 
Metropolis scheme, the acceptance of such trial moves depends on the ratio of the Boltzmann 
weights of the new and old states for each temperature T lfT8ll . However, if this were the only 
criterion, there would be a tendency to generate homopolymer chains with a low energy, rather 
than chains that fold selectively into a specific target structure. To ensure a type composition far 
from the homopolymer region of sequence space, we sample the equilibrium sequences for the 
generalized energy function 

W = E-z p lnNp (5) 

where E is the energy of the patchy polymer defined as the sum of all the contributions calculated 
according to Eq.|2](as the structure does not change the spring and directional terms are constant), 
£ p is a scale factor adjusting the relative importance of the two terms in the equation, and Np is the 
number of permutations that are possible for a given set of particle of given types, 

N\ 

Np= . . . (6) 

n\ \n2\ny. - -hmI 

Here, N is the total number of particles in the chain, M is the number of different particles types, 
and wi ,n2, • • • -> n M, etc. are the numbers of particles of type 1 , 2, • • • ,M, respectively. In the present 
study we have considered an alphabet of only two particles types, namely H and P. Hence, it is 
simply the ratio: 

N\ 

N P =——. (7) 
hh \np\ 
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While sampling sequence space with the Monte Carlo scheme, we set £ p to a sufficiently high 
value, £ p — 5k^T, to generate sequences with a heterogeneous composition. Note that the ad- 
ditional term to the generalized energy of Eq. ([5]) is a phenomenological bias that drives the 
simulation towards the region of heterogeneous sequences. Although due to this particular bias 
we may miss some low energy configurations, we believe that an exhaustive sampling of the entire 
sequence space of each structure is not necessary to compare the designability. 

During each simulation we projected the configurational space on the collective variable Qs, 
which counts the number of contacts between the spheres. We then compute the free energy F(Qs) 
as a function of Qs, 

F(Qs) = -kTln[P(Q s )]. (8) 

Here, P(Qs) denotes the normalized histogram of the number of sampled conformations with 
order parameter Qs. In practice, a direct calculation of this histogram is not efficient, since even 
such short chains tends to be trapped in local minima, especially at low temperatures. To induce 
escape from these local minima, we made use of the Virtual Move Parallel Tempering Monte 
Carlo sampling scheme proposed by Coluzza and Frenkel lfT9L based on the Waste Recycling 
approach EDll . This scheme is very efficient in sampling both high and low free energy states. 



III. RESULTS 



We first carried out a simulation in which we varied both the configuration as well as the se- 
quence of the chain for the patch-patch interaction PI from Eq.[3} We considered the combination 
of parameters in Table [I] and we projected the structure space on the collective variable Qs, which 
is a count of the number of contacts between the spheres. In Fig. [2j we show the free energy 
profiles F(Qs)/kftT as a function of Qs at various combinations of the potentials ranges r max and 
Gpi. From the plot we can see that for r max = 2.2, 2.5 and Op\ = 0.2 the landscapes present a sin- 
gle sharp minimum, indicating that there is a preferred number of contacts for the structures that 
have high occurrence. As r max is increased, the position of the minimum, not surprisingly, shifts 
towards higher values of Qs and it also widens suggesting a richer structural space. We further 
characterized the ensemble of configurations by looking at the structures that correspond to the 
global minimum of each free energy profile. In Figs. [3] and [4] we show a real space representation 
of such structures. The first aspect that becomes apparent is that there is a wide range of possible 
structures (it is important to remember that the structures in the minimum are only the most com- 
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Figure 2: Free energy F{Qs)/k&T as a function of the number Q$ of contacts among the spheres obtained 
for a patchy polymer of length N = 20 with 3 patches per particle interacting with the directional interaction 
defined in Eq. |3j The simulations were performed by altering the structure of the chain as well as by 
mutating the identities of the spheres that intervene in Eab interaction in Eq.[2| at the temperature T = 0.2. 
Each curve represents a simulation with a different pair of values (r max ,0/>i) as indicated in Table [i] To 
each pair corresponds a different symbol, however we have grouped the curves with the same range of the 
directional potential E P \ defined in Eq. [3] 

mon and not the only possible ones). In addition, we observed the appearance of double helices 
at the values r max and Cp\ for which the interacting patches sit very close to the half distance r max 
between the centers of the spheres ( the case r max = 2.2, Gp\ = 0.2 is a bit of an exception as it 
shows only a partial helical configuration). It is important to notice that the pitch of the double 
helical structures in Fig. [4] shrinks with increasing values of r max . One could imagine a system 
where the range of the isotropic interaction is controllable by external parameters (e.g. PH, DNA 



9 



linkers, depletion), which in turns will allow for the control over the extension of the helix. 




Figure 3: Typical configurations obtained from the simulations of the interaction potential PI (Eq.[3]). For 
clarity we have not represented each sphere along the chain, but instead we draw only the backbone as a thick 
line joining the centers of the spheres. The patches are represented as white stick anchored on the center of 
the colloid and of length equal to the radius of the colloid, so that the end of the stick corresponds to the 
position in space of the patch. The structures have been taken from the ensemble of structures corresponding 
to the minimum of the free energy F{Qs)/k^T for the combinations of parameters: a) r max = 2.2 o P \ = 0.2, 
b) r max = 3.0 an = 0.2, c) r max = 4.0 o P1 = 0.2, d) r max = 6.0 a P1 = 0.2, e) r max = 6.0 o P1 = 1.0. The 
configurations have a loop structure that folds on itself very similar to what in proteins is called a beta- 
sheet. We colored structure (a) in red to highlight the loop-helix nature of the backbone which resembles 
the helices in Fig. [^corresponding to the other set of parameters. 

We now consider the second interaction potential P2 from Eq. [4} We repeated the same pro- 
cedure described above for the PI interaction potential. In Fig. [5} we plot the free energies 
F{Qs)/k%T for the values of the parameters r max , G P 2 in Table [l| as a function of the same collec- 
tive variable Qs. As before we observe a clear trend that the value of Qs at the minimum increases 
for increasing values of r max . However, if we look at the actual structures we noticed that they are 
all very close to each other (see Fig. [6]) indicating that the displacement of the minimum is just an 
effect of counting more neighbors per particle. This is true also for the longer ranges of the direc- 
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Figure 4: As for the calculations shown in Fig. [3] we have isolated the characteristic structures of the 
minimum free energy F(Qs)/k&T but now for the combinations of parameters: a) r max = 2.5 Gp\ = 0.2, 
b) r max = 4.0 Gp\ = 1.0, c) r max = 6.0 Gp\ = 2.0. All the systems exhibit a strong preference for helical 
configurations of the backbone. To be noted is the shrinking effect that helices show upon increase of the 
range r max of the isotropic interaction. 

tional potential, that show for Op2 = 3.0 a double minimum and for Gp2 = 4.0 a very wide basin. 
The latter is the result of large fluctuations that the long range potentials can accommodate. It is 
striking to obtain such a dramatic preference for a specific configuration just by slightly altering 
the form of the interaction potential. By using the same form of the potentials for both isotropic 
and directional contributions we have reduced the frustration between the two potentials, allowing 
the system to find an optimal structure for all sets of parameters. We stress that we have chosen 
the values of the parameters in order to fix the second viral coefficient of all the terms of the po- 
tential constant across all simulations. Slightly different configurations occur when the range of 
the directional interaction Gp2 is half of range r max the isotropic interaction (See (e) and (g) of 
Fig. [6]). The main difference between the configurations (e) and (g) and the others is in the central 
loop of the backbone, that for larger values of Gp2 expands into an extended arch. In particular the 
arch allows for large compressions of the molecule that results in the wide minima for Gp2 = 3.0 
in Fig. [5} 
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Figure 5: Free energy F{Qs)/k^T as a function of Qs for the same simulation conditions and parameters as 
in Fig.[2j but with the patch-patch interaction from Eq.[4j Each curve represents a simulation with a different 
pair of values (r max ,0/>2) as indicated in Table [i] To each pair corresponds a different symbol (see inset), 
however we have grouped the curves with the same range of the directional potential E P 2 defined in Eq. |4j 
The wide free energy minima observed for the r max = 6.0, G/>2 = 4.0 case is due to the large fluctuations 
that the long range potentials can tolerate. 

IV. SUMMARY AND DISCUSSION 

In this work, we present a methodology to explore the configurational space of a chain of col- 
loidal particles decorated with a single patch. The particles interact through an isotropic potential 
that depends on the identity of the particles as well as through a direction dependent patch-patch 
interaction. We developed a methodology to sample the probability of observing a structure with 
many different sequences therefore allowing for a comparison of the relative probability of observ- 
ing each configuration. The complexity of such a space gives an idea of how a particular model is 
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Figure 6: Minimum free energy configurations (see Fig. [5]), obtained from the simulations of the P2 (Eq.[4]) 
interaction potential. Each structure is representative of different combinations of the potentials range: 
a) r max 2.2 a P2 2.2, b) r max 2.5 o P2 2.2, c) r max 3.0 o P2 2.2, d) r max 4.0 o P2 2.2, e) 
^max = 4.0 Gp2 = 2.0, f) r max = 6.0 Gp2 = 2.2, g) r max = 6.0 <T/> 2 = 3.0. We colored structures (e) and (g) 
in blue because they correspond to the same set of parameters that for the PI interaction potential gave rise 
to the structures in Fig. [4] Unexpectedly all the structures are almost identical (with an average root mean 
square distance between the centers of the spheres of 1 Rhc)- The common backbone arrangement appears 
to be a distorted helix where the patches interact via the tails of the potential in Eq. [4] 

suitable for the design of self-assembling structures. We then defined two different forms for the 
interaction potential and for each of the potentials we considered a wide range of values for the 
parameters of the potential. From the analysis of the projection of such a space on two collective 
variables, we identified a wide range of structures with high probability of being observed. 

Our analysis demonstrates that the combination of interactions, with which the chain of col- 
loids is decorated, is capable of shaping the configurational space to a reduced ensemble of con- 
figurations. Furthermore, we have shown that the reduced ensemble can be dramatically altered 
by changing few parameters of the interaction potential (e.g., the interaction range). Such con- 
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trol over the size and the structure heterogeneity of the configurational space makes the chain of 
patchy-colloids an ideal candidate for the experimental realization of new self-assembling systems. 
Moreover, we have provided initial but strong indications that, by inducing frustration between the 
directional (Eq. [3]) and the isotropic potential the chain is forced to adopt a wider range of struc- 
tures, compared to the scenario where both potentials shared the same functional form (Eq.[4]). In 
particular, for specific parameters the backbone of the chain has a strong preference for ordered 
double helices with a pitch controllable by the range of the isotropic interaction. 

Our results provide an important starting point for the realization of designable chains of patchy 
colloids. The methodology used here can be easily extended to an arbitrary set of interactions and 
provide, as a result, a set of target structure for the design of chains of patchy colloids. 
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